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Memory effects play a key role in the dynamics of strongly correlated systems driven out of equi- 
librium. In the present study, we explore the nature of memory in the nonequilibrium Anderson 
impurity model. The Nakajima-Zwanzig-Mori formalism is used to derive an exact generalized 
quantum master equation for the reduced density matrix of the interacting quantum dot, which 
includes a non-Markovian memory kernel. A real-time path integral formulation is developed, in 
which all diagrams are stochastically sampled in order to numerically evaluate the memory kernel. 
We explore the effects of temperature down to the Kondo regime, as well as the role of source-drain 
bias voltage and band width on the memory. Typically, the memory decays on timescales signifi- 
cantly shorter than the dynamics of the reduced density matrix itself, yet under certain conditions 
it develops a smaller long tail. In addition we address the conditions required for the existence, 
uniqueness and stability of a steady-state. 



Interest in the problem of intrinsically nonequilibrium 
open quantum systems — in which one considers a small, 
strongly interacting and highly correlated region coupled 
to several large, noninteracting baths — has been surging 
in both experiment and theory. The aim of theory in 
this regard is to provide a solid framework for under- 
standing phenomena ranging from the nonequilibrium 
Kondo effect in quantum dots to conductance through 
single molecules [T]. While much progress has been 
made recently, based on brute-force approaches such as 
time-dependent numerical renormalization group tech- 
niques [2H1] and iterative |5J |5] or stochastic [7HH] dia- 
grammatic approaches to real-time path integral formu- 
lations, the problem has never been fully solved in a sat- 
isfactory manner. In fact, it is becoming clearer that ma- 
jor gaps exist in our understanding of the dynamics, the 
crossover regimes, the dependence on initial conditions 
and the behavior at steady state. 

The kind of open systems discussed above are often 
addressed by impurity models, which explicitly account 
for the two types of regions within the problem by par- 
titioning the Hamiltonian into system and bath sub- 
spaces: H = Hs + Hb + V, where H$ represents a low- 
dimensional but interacting "system" subspace, Hb rep- 
resents a set of non-interacting lead or bath subspaces, 
and V is a system-bath coupling term. The dynamics 
generated by such Hamiltonians can feature transients on 
timescales that are much longer than the typical inverse 
energy scale |10| , where numerically exact approaches be- 
come intractable due to the exponential growth of the ac- 
tive space or the equivalent complications resulting from 
the dynamical sign problem. In many important situa- 
tions, however, the non-interacting baths can be traced 
out, leading to a reduced description of the dynamics of 
the interacting system at the cost of introducing non- 
locality in the time propagation In path integral 
approaches the effects of the leads are accounted for by 
a time non-local influence functional [5) [7] . 



Perhaps a more appealing approach, which has been 
used to derive very successful perturbative schemes for 
fermionic systems |121 fT5] but is notoriously difficult to 
carry out exactly, is based on the generalized quantum 
master equation (GQME). In this formalism, a so-called 
memory kernel replaces the influence functional. The 
complexity of solving the many-body quantum Liouville 
equation is then reduced to the evaluation of this memory 
kernel, which fully determines the dynamics of the sys- 
tem. Furthermore, the memory kernel contains all the 
information needed to resolve questions concerning the 
existence and uniqueness of a steady-state |141 [T5] . as 
well as the values of system observables at steady state. 
While the dynamical timescale of the system typically ex- 
ceeds the characteristic inverse energy scale, the memory 
kernel is expected to decay on relatively short timescales 
for a large and interesting class of physical situations (es- 
sentially whenever the bandwidth of the baths is much 
larger than other energy scales in the problem). Thus, 
brute-force approaches limited to short times are well 
suited to its calculation. Once the memory function is 
known, the formalism is exact and tractable at all times. 

In this Letter, we explore the nature of memory in 
nonequilibrium impurity models, and focus on the An- 
derson problem [H], covering the effects of temperature 
down to the Kondo regime, the role of source-drain bias 
voltages and band width. This is accomplished by adopt- 
ing the Nakajima-Zwanzig-Mori |17ffT§] formalism to 
derive an exact GQME for the reduced density matrix 
a (t) = Pp (t) of the interacting system, which includes a 
non-Markovian memory kernel. The conjecture that the 
memory decays on timescales significantly shorter than 
the dynamics of a (t) is confirmed, yet it is found that 
it develops a smaller long tail when the Hubbard term 
is switched on. The approach provides means to simu- 
late the dynamics of the strongly correlated subsystem on 
timescales beyond the limits of the path integral method 
itself and reveals the conditions required for the exis- 
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Figure 1. Distinct nonzero memory kernel elements for an 
initially unoccupied dot and several values of the interaction 
energy U. The remaining parameters are 4r- = — = s, 
f? = 20, and /3T = 1. 

tence, uniqueness and stability of a steady-state under a 
finite source-drain bias. 

The exact equation of motion of the complete density 
matrix, ih^p — [H, p] , is governed by the Liouvillian L = 
[H, ...). If the coupling term V were to be turned off, the 
dynamics of the two subspaces would be given by similar 
equations with the Liouvillian replaced by the system 
and bath Liouvillians, £$ — [Hs, ...] and C B = [H B , •■•], 
respectively. We also define a coupling Liouvillian Cv = 
[V, ...]. The equation of motion for the reduced density 
operator pga = Pp (using the projection operator onto 
the system subspace P = pbTyb {•••}) is then given by: 

ih^ t a{t)^L s o{t)+d{t)- l -J^ dT K (r)a(t-T), (1) 
where (i) = Tr B {c v e-i Q£t Qf>o} 

contains the initial 

correlations, k (t) = Tib {£ve~ i ® Cl ~Q£-Pb} is the mem- 
ory kernel in superoperator form, and Q = 1 — P. The 
initial condition are contained within the initial density 
matrix po, which also determines the initial bath part of 
the density matrix ps = Qpo- 

The above equation for a (t) is exact, yet requires as in- 
put the two superoperator functions k (t) and i9 (t), both 
of which have been defined in terms of projected propa- 
gation. Evaluating such projected dynamics is cumber- 
some, and can be reduced to solving a superoperator 
Volterra equation of the second type involving a quan- 
tity $ (t) which is free of projected propagation |2"Ul I21j . 



For the memory kernel one finds: 

K (t) = ih$ (t) - $ (r) C s 

+±J dr$ (*-t)k(t), (2) 

^(t) = TT B {c v e-i ct p B }. (3) 

A similar procedure exists for the initial correlation term 
however, here we consider only the initially factorized 
case po = pb <£> c(0), for which d — 0. As discussed 
below, the matrix elements of the superoperator <!> are 
identical to quantities to which RT-PIMC has already 
been applied [7J |S] . 

In the Anderson impurity model Hs — 2^j = +i £id\di + 
Ud\d t d\d l7 H B = T,k,i=U £ika ik a ik and v = 

E/si=t4.*»fe^* a lfe + **k°*k<4- Thus, the system subspace 
is four-dimensional, being spanned by the states = 
1 00) = |0) , 1 = |01) = 4 1°) » 2 = l 10 > = 4 1°) and 
3 = 1 11) = d\d\ |0). With this notation, we can per- 
form a calculation to derive an expression for the system 
Liouvillian 

+£2 {Sqi + 5 q3 - 5ji - 8js) 
+U(S q3 -S j3 )}, (4) 

and for which takes the more complicated form 

<p$ =2&^T* B {p B tf\AF>\q)}, (6) 
k 

^ = -2j2Tr B {p B ( q '\BW\ q )}, (7) 

k 

where A k x) = ti k did 2 d\a\ k , A k 2) = h k did\d 2 a\ k , A k 3) = 
t 2 kd\d\d 2 a\ kl a£> = t 2 kd\did 2 a\ k , B k x) = t lk d 2 a\ k and 
B£' = t 2k d\a} 2k . All the quantities in ([5])-([7]) are implic- 
itly time-dependent and can be evaluated directly with 
RT-PIMC 

One can draw a few analytical conclusions directly 
from the block structure of First, while both the A k m ^ 
and B k m ^ operators introduced above conserve the to- 
tal particle number, only the A k conserve the particle 
number for each spin. Therefore, the are nonzero 

only if q,q' = 1,2 or 2,1 (since states 1 and 2 are the 
only dot states which have the same total occupation, 
but differ in per-spin occupation), while the <p£j can be 
nonzero only when q = q'. From ^ we can then imme- 
diately see that the diagonal density matrix elements arc 
coupled only to each other, with the two singly occupied 
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Figure 2. Total population derivative from direct RT-PIMC 
data compared with the results of the memory-kernel formal- 
ism (left panels) and predicted dot populations (right panels) 
for an initially unoccupied dot for several values of the inter- 
action energy U . Remaining parameters are the same as in 
Fig.jl] The cutoff time is shown on the right panels as a 
vertical dotted line. 



off-diagonal elements |1) (2| and |2) (1| forming a second 
closed block. If our interest is limited to the state pop- 
ulations only the 16 (instead of 64) functions <fqq need 
be evaluated within RT-PIMC. 

In Fig. [I] we plot the nonzero elements of the mem- 
ory kernel for different values of U. Due to the 
block structure and the symmetric choice of parame- 
ters, only seven distinct nonzero elements exist (two 
for {7 = 0). To make the parametrization definite 
within the simulations, we assume lead coupling den- 
sities of the form T. l (E) = 2nJ2 k \tk\ 2 S{E~e lk ) = 

( 1+e ,( E -ec))(i +e -.( g+ e e )) ' where £ c is the band cutoff en- 
ergy and v is the inverse of the cutoff width. T is the 
maximum value attainable by T (E) — J^. (E), as well 
as its wide band limit, and will be the unit of energy 
throughout the following text. We also concentrate on 
£j = — 7^, known as the symmetric case of the Anderson 
model, yet the formalism is general and this is certainly 
not a requirement The temperature (3 and the chemical 
potentials /il and /i^ enter the calculation through pg at 
time zero for which we assume the proper grand canonical 
distribution. 

In the noninteracting case (top panel of Fig. [I]), a rapid 
decay to zero is observed. Despite the relatively broad 
and soft-edged band chosen here, the decay occurs over 
a timescale smaller than the inverse coupling but compa- 
rable to it. We can learn from this that any approxima- 
tion based on short memory should be expected to fail 
unless it can allow for memories of at least this length, 
meaning for instance that Markovian approaches to the 



problem cannot be expected to succeed in general. As 
U is increased, it becomes clear that the interaction — 
despite breaking most of the symmetries between the 
various elements — does not significantly affect memory 
decay on the first timescale. However, even at very small 
interaction energies, a second, longer timescale develops: 
at this timescale, a small part of each memory kernel 
element decays more slowly. 

The formalism becomes extremely interesting if hav- 
ing the memory as an input only up to some finite cut- 
off time t c — at which the system's dynamics have not 
yet died out — allows accurate predictions at far longer 
times. This will occur if the memory function has essen- 
tially gone to zero by this time, such that it can be safely 
truncated. In the left panels of Fig. [2] we plot the time 
derivative of the total population (dP/dt) and show that 
for certain noninteracting parameters this does indeed 
happen: once Tt c > | dynamics at times over an order 
of magnitude greater than those of the memory kernel 
are reliably reproduced, and the exact steady-state re- 
sult for all diagonal elements of a is obtained to within 
the numerical errors and shown on the top right panel 
of Fig. [2] However, as might be expected from the anal- 
ysis of Fig. [T] in the strongly correlated cases a short 
cutoff time results in incorrect populations when prop- 
agated for much longer times than t c , since truncating 
the memory at that point has not yet become physically 
reasonable. For the cases of U = T and U = 6T the qual- 
itative physics is captured correctly within t c = 3/2T 
in that depopulation of the zero-electron level is acceler- 
ated at short times; for U — 6T one also observes that at 
longer times the one-electron levels draw most of the pop- 
ulation while the more energetic zero- and two-electron 
levels are suppressed. However, for U = 3r, the behavior 
predicted by the truncated memory function is clearly 
wrong, consistent with the existence of long time tails in 
the memory. 

The top two panels of Fig. [3] explore the decay of the 
memory kernel elements more clearly by plotting the av- 
erage absolute value of the memory kernel elements on 
a logarithmic scale, for various values of U, e c and j3. 
Notably, while only the noninteracting case appears to 
have a memory kernel which goes to within the numeri- 
cal errors of zero within the simulation timescale shown, 
strong interaction actually appears to reduce the memory 
lifetime when compared with intermediate values. This 
relates to the fact the U — ST problem is "harder" in 
this sense than the strongly interacting U = 6T problem, 
as discussed above. While increasing either the band 
width or the temperature appears to affect the short- 
term behavior, reducing the shorter memory timescale, 
the longer timescale appears largely unaffected by the 
variation of these parameters. This observation seems to 
be consistent with the hypothesis that the timescale of 
the tail's memory decay is related to the inverse Kondo 
temperature. However, we find that a large bias voltage 
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Figure 3. The average absolute value of memory kernel el- 
ements at V± = = |, ^ = 20, /?r = 1 and different 
values of the interaction energy U (top panel); the same at 
fJ-L = Mfl = 0, y~ =40, y ~6 and different temperatures /3 
(middle panel); and the predicted steady state values of the 
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does not affect markedly the longer timescale, despite 
supposedly destroying the Kondo correlations. 

In addition to the time dependence, one can obtain 
the steady state result directly from the stationary state 
equation 



Cs — —k(z—>i0) 



oo 



lim z-d (z) (8) 

2— >0 



with the added condition that Trcr = 1. This means 
that for an initially uncorrelated system, a unique 
steady state exists if and only if the supermatrix 

Tr s [C s - \h{z -> »0)] \k) has a degener- 

acy of exactly one. In the bottom panel of Fig. [3] the 
steady state values obtained from this formula for pa- 
rameters close to the Kondo regime are plotted against 
the cutoff time t c . While the trace of the density matrix is 
conserved, physically impossible results appear at inter- 
mediate cutoff times and convergence is not yet achieved. 
The long tails of the memory kernel elements are there- 
fore crucially important for the correct prediction of both 
dynamical and steady-state properties. 

In conclusion, we have developed a numerically ex- 
act method for the formulation and solution of the re- 
duced dynamics of quantum impurity models, and ap- 
plied it to the nonequilibrium Anderson model. It is clear 
from our results that the physics of even a noninteract- 



ing electronic junction cannot be fully captured within 
a Markovian picture, and that on-site interaction results 
in deeply non-Markovian physics even at relatively large 
band widths, bias voltages and temperatures. We show 
that the long memory tails induced by the Hubbard term 
affect both the dynamics and the steady-state, despite 
their relative smallness. In the computational sense, the 
proposed method is extremely useful in extending the 
applicability of RT-PIMC to long time scales and steady 
state when the memory goes to zero within the simula- 
tion time scale, but the dynamics do not. 
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